Modeling the Two-Body Problem (Earth-Sun)¶

# Two-Body Problem (Earth–Sun) — Python Notebook Template (ELLIPSE VERSION)

# -----------------------------------------------------------------------

# Change: set v0 to be a small factor times the circular speed instead of a fixed number.

​

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

​

# ----------------------------

# 1) Constants (SI units)

# ----------------------------

G = 6.67430e-11  # m^3 / (kg s^2)

​

m_sun = 1.98847e30     # kg

m_earth = 5.9722e24    # kg

​

mu = G * (m_sun + m_earth)  # gravitational parameter for relative motion

​

# ----------------------------

# 2) Initial conditions

# ----------------------------

r0 = np.array([1.496e11, 0.0, 0.0])     # meters (about 1 AU)

​

# --- CHANGED: compute circular speed and scale it ---

r0_mag = np.linalg.norm(r0)

v_circ = np.sqrt(mu / r0_mag)           # circular speed at r0

v_esc  = np.sqrt(2 * mu / r0_mag)       # escape speed at r0 (for reference)

​

factor = 0.85  # <--- CHANGED: try 0.995, 0.99, 1.005, 1.01

v0 = np.array([0.0, factor * v_circ, 0.0])

​

print("v_circ =", v_circ, "m/s")

print("v0     =", np.linalg.norm(v0), "m/s  (factor =", factor, ")")

print("v_esc  =", v_esc, "m/s")

​

y0 = np.hstack([r0, v0])  # state = [r, v] in R^6

​

# ----------------------------

# 3) ODE: relative two-body

# ----------------------------

def two_body_relative(t: float, y: np.ndarray, mu: float) -> np.ndarray:

    r = y[:3]

    v = y[3:]

​

    r_norm = np.linalg.norm(r)

    if r_norm == 0:

        raise ValueError("Encountered r_norm = 0 (collision / singularity).")

​

    a = -mu * r / r_norm**3

    return np.hstack([v, a])

​

# ----------------------------

# 4) Simulation horizon

# ----------------------------

years = 2

​

T = years * 365.25 * 24 * 3600  # seconds

t_span = (0.0, T)

​

n_points = 5000

t_eval = np.linspace(t_span[0], t_span[1], n_points)

​

sol = solve_ivp(

    fun=lambda t, y: two_body_relative(t, y, mu),

    t_span=t_span,

    y0=y0,

    method="DOP853",

    rtol=1e-10,

    atol=1e-12,

    t_eval=t_eval,

)

​

if not sol.success:

    raise RuntimeError(f"ODE solver failed: {sol.message}")

​

t = sol.t

Y = sol.y.T

​

r = Y[:, :3]

v = Y[:, 3:]

​

# ----------------------------

# 5) Plots: orbit + radius

# ----------------------------

plt.figure()

plt.plot(r[:, 0], r[:, 1])

plt.gca().set_aspect("equal", adjustable="box")

plt.grid(True)

plt.xlabel("x (m)")

plt.ylabel("y (m)")

plt.title(f"Two-body relative orbit (Earth–Sun), ellipse (factor={factor}), {years} year(s)")

plt.show()

​

​

# --- pick one time index (0 is fine) ---

i = 0

r0_vec = r[i]

v0_vec = v[i]

​

r0_mag = np.linalg.norm(r0_vec)

v0_mag = np.linalg.norm(v0_vec)

​

# Specific angular momentum vector

h0 = np.cross(r0_vec, v0_vec)

​

# Eccentricity vector: e = (v x h)/mu - r/|r|

e_vec = (np.cross(v0_vec, h0) / mu) - (r0_vec / r0_mag)

e = np.linalg.norm(e_vec)

​

print("eccentricity e =", e)

print("e_vec =", e_vec)

​

# --- Build a rotation matrix so that e_vec points along +x ---

ex = e_vec / np.linalg.norm(e_vec)              # unit vector along eccentricity direction

hz = h0 / np.linalg.norm(h0)                    # unit vector normal to orbital plane

ey = np.cross(hz, ex)                           # completes right-handed basis

​

# Rotation matrix: columns are new basis vectors in old coordinates

R = np.vstack([ex, ey, hz])  # (3x3) maps old -> new via R @ r_old

​

# Rotate all positions into the ellipse-aligned frame

r_rot = (R @ r.T).T  # shape (N,3)

​

plt.figure()

plt.plot(r_rot[:,0], r_rot[:,1])

plt.gca().set_aspect("equal", adjustable="box")

plt.grid(True)

plt.xlabel("x' (m) (along eccentricity vector)")

plt.ylabel("y' (m)")

plt.title("Orbit in ellipse-aligned frame (ellipse should be obvious)")

plt.show()

​

​

​

​

r_norm = np.linalg.norm(r, axis=1)

​

plt.figure()

plt.plot(t, r_norm)

plt.grid(True)

plt.xlabel("t (s)")

plt.ylabel("||r(t)|| (m)")

plt.title("Distance over time (ellipse shows clear min/max)")

plt.show()

​

print("r_min =", r_norm.min())

print("r_max =", r_norm.max())

print("approx eccentricity from min/max:",

      (r_norm.max() - r_norm.min()) / (r_norm.max() + r_norm.min()))

​

# ----------------------------

# 6) Long-run behavior checks

# ----------------------------

v_norm = np.linalg.norm(v, axis=1)

​

energy = 0.5 * v_norm**2 - mu / r_norm

​

h = np.cross(r, v)

h_norm = np.linalg.norm(h, axis=1)

​

plt.figure()

rel_energy_drift = (energy - energy[0]) / abs(energy[0])

plt.plot(t, rel_energy_drift)

plt.grid(True)

plt.xlabel("t (s)")

plt.ylabel("relative energy drift")

plt.title("Relative energy drift")

plt.show()

​

plt.figure()

rel_h_drift = (h_norm - h_norm[0]) / abs(h_norm[0])

plt.plot(t, rel_h_drift)

plt.grid(True)

plt.xlabel("t (s)")

plt.ylabel("relative |h| drift")

plt.title("Relative angular momentum drift")

plt.show()

​

​

​

from matplotlib.animation import FuncAnimation

from IPython.display import HTML

​

def animate_distance_plot(t, r_norm, step=5):

    """

    Animates distance vs time using FuncAnimation.

    step controls how many points are added per frame.

    """

    t_years = t / (365.25*24*3600)

​

    fig, ax = plt.subplots()

    ax.set_xlabel("t (years)")

    ax.set_ylabel("||r(t)|| (m)")

    ax.grid(True)

​

    line, = ax.plot([], [], lw=2)

​

    # Set axes limits once

    ax.set_xlim(t_years[0], t_years[-1])

    ax.set_ylim(r_norm.min(), r_norm.max())

​

    def init():

        line.set_data([], [])

        return (line,)

​

    def update(frame):

        k = frame * step

        if k < 2:

            k = 2

        x = t_years[:k]

        y = r_norm[:k]

        line.set_data(x, y)

        return (line,)

​

    n_frames = int(np.ceil(len(t) / step))

    anim = FuncAnimation(fig, update, frames=n_frames, init_func=init, blit=True, interval=30)

​

    plt.close(fig)

    return HTML(anim.to_jshtml())

​

# Example usage:

animate_distance_plot(t, r_norm, step=10)

​